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Cardiovascular networks span the body by branching across many generations of vessels. The 
resulting structure delivers blood over long distances to supply all cells with oxygen via the relatively 
short-range process of diffusion at the capillary level. The structural features of the network that 
accomplish this density and ubiquity of capillaries are often called space-filling. There are multiple 
strategies to fill a space, but some strategies do not lead to biologically adaptive structures by 
requiring too much construction material or space, delivering resources too slowly, or using too 
much power to move blood through the system. We empirically measure the structure of real 
networks (18 humans and 1 mouse) and compare these observations with predictions of model 
networks that are space-filling and constrained by a few guiding biological principles. We devise 
a numerical method that enables the investigation of space-filling strategies and determination of 
which biological principles influence network structure. Optimization for only a single principle 
creates unrealistic networks that represent an extreme limit of the possible structures that could be 
observed in nature. We first study these extreme limits for two competing principles, minimal total 
material and minimal path lengths. We combine these two principles and enforce various thresholds 
for balance in the network hierarchy, which provides a novel approach that highlights the trade-offs 
faced by biological networks and yields predictions that better match our empirical data. 

PACS numbers: 87.10.Vg 89.75.Da 89.75.Fb 89.75.Hc 89.75.Kd 


I. INTRODUCTION 

The vital functions of the cardiovascular system are the 
distribution of oxygen and nutrient resources throughout 
the body, as well as the collection and filtration of waste 
by circulating blood. Transfer of resources and waste oc¬ 
curs primarily at the capillary level via diffusion through 
nearby tissue. Consequently, this smallest level of the 
network must reach all living cells in order to maintain 
them, filling the entire space of the body. In models 
developed by Krogh for effective diffusion of oxygen [lj], 
cells cannot survive beyond a maximum distance from a 
capillary. This defines a service volume of cells that is as¬ 
sociated with each capillary, which has a typical size that 
has been observed to vary across species based on cellular 
metabolic rate [2, 3]. The constraint on maximum dis¬ 
tance from capillaries necessitates that the final levels of 
the cardiovascular network are also space-filling through¬ 
out the body. In this paper we investigate the relation 
between this space-filling property and basic optimiza¬ 
tion principles such as the minimization of costs from 
construction material and pumping power. Specifically, 
we highlight how this relation influences the asymmetries 
in sizes and flow rates of sibling segments as measured in 
empirical data. 

A central focus of our investigation of cardiovascular 
systems is the space-filling properties of networks, but 
these properties are also of great interest in many other 
contexts. General space-covering hexagonal patterns ap¬ 
pear in nature in the cell structure of beehives as well 
as in economic theories for market areas Q. Trees (the 
woody, perennial plants) have been studied for both how 
forests fill an area [Bj], as well as how the vascular struc¬ 
ture within an individual plant determines the hydraulics 


of resource delivery §3, 7]. Apollonian networks [8] de¬ 
scribe the space-filling packing of spheres of various sizes, 
which are similar in the cardiovascular system to consid¬ 
ering the volumes of different subtrees of the network. 
Efficiently filling space in two dimensions is important 
for information visualization [9]. In addition to these ap¬ 
plications, Kuffner and LaValle study space-filling tree 
networks (i.e., networks that branch and have no closed 
loops) to determine a route from one location to another 
jlQj |. but without the biological constraints that we im¬ 
pose here. For cardiovascular networks, this motion plan¬ 
ning is analogous to how the vascular structure routes 
blood. Efficient routing of blood through a hierarchy is 
central to models that investigate allometric scaling of 
metabolic rate with body mass [2|, 3, 11-[L4[, which build 
on metabolic scaling by Kleiber |15|. Determining spe¬ 
cific space-filling strategies will inform these models to 
better describe the cardiovascular system. 

Developmental processes (i.e., growth) as well as evo¬ 
lutionary pressures, such as efficiency in material and 
energy use, shape the structure of cardiovascular net¬ 
works. Filling a volume or surface efficiently with the 
terminal nodes of a branching network is nontrivial, es¬ 
pecially when the distribution system must reliably de¬ 
liver blood at each stage of development. The system 
must also be robust to changes in vessel lengths and the 
number of hierarchical levels that result from growth and 
obstructions from damage US S3- We propound that 
these structural challenges lead to the asymmetric fea¬ 
tures in both the local relations of segments, as well as 
in the global paths from the heart to each service volume 
that we observe in empirical data. 

At the local level, the ratio of lengths between par¬ 
ent and child segments may vary across the network, de¬ 
viating from strict self-similarity. Additionally, sibling 
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segments may vary in length, which we call asymmetric 
branching. Within our numerical method, these features 
arise as a result of optimizing branch point positions rel¬ 
ative to adjacent branch points across the network. Vari¬ 
ation in the length and number of levels in paths means 
that the network is also not symmetric or balanced in 
these global properties. By allowing these asymmetries 
and explicitly ensuring space-filling structure, we expand 
other models that are strictly balanced in network hier¬ 
archy and perfectly symmetrical in local quantities. 

Asymmetries observed in real systems motivate our 
investigation of the space-filling properties and asym¬ 
metries in cardiovascular networks. These observations 
show that such networks have a tendency for the lengths 
of sibling segments to be distributed less symmetrically 
than is the case for radii. The empirical data in Sec. [V| 
come from two types of sources. Images collected through 
microtomography of mouse lung comprise one data set. 
The mice were part of a study with different manipula¬ 
tions of matrix GLA protein (which causes the vascula¬ 
ture to be under- or over-branched 18]), but we focus 
on the wild-type specimen for our analysis. The other 
data set, collected through MRI, excludes the lungs and 
instead focuses on the central vasculature in the human 
head and torso [19]. We utilize the custom software Angi- 
cart M to analyze these two distinct vascular data sets. 
Because of the noninvasive nature of these data acqui¬ 
sition and analysis techniques, future studies have the 
opportunity to track the development of cardiovascular 
systems as individual organisms grow and age, including 
repair after the system incurs damage (i.e., wound heal¬ 
ing). Such studies can elucidate the effects that patterns 
of growth and changes from damage have on the final, 
mature state of the network. 

In this paper, we study the optimization principles that 
correspond to evolutionary pressures for efficiency in ma¬ 
terial cost and power loss. Our focus is the influence of 
space-filling patterns on length asymmetry distributions 
without the explicit inclusion of radius information. The 
list of candidate networks includes all possible hierar¬ 
chical (topological) connections between the heart and 
all capillaries. For each hierarchy and unique permu¬ 
tation of pairings between terminal vessels and service 
volumes (see Fig. fid]) , we must determine the positions 
of the branch points. The combination of the hierarchy, 
service-volume pairings and branch point positions de¬ 
fines the configuration of a candidate network. For these 
reasons, we must search through many candidate config¬ 
urations to determine the most efficient structure. We 
quantify the fitness of each candidate network using in¬ 
dividual segment lengths between branch points as well 
as full path lengths between each capillary and the heart. 

To perform a reliable comparison between candidate 
configurations, it is crucial to determine branch point 
positions in a consistent way. We determine these posi¬ 
tions iteratively for the entire network in order to identify 
the global optimum. While the local process of choosing 
branch points that minimize total vessel lengths (or sim¬ 


ilar features) is relatively straightforward to iterate over 
the network, any single branch point and its relation to 
its neighbors relies indirectly on updates that are applied 
elsewhere. This dependence emerges from the fact that 
each end of a vessel is connected to a branch point, which 
upon repositioning affects the lengths of all vessels that it 
joins. The uniqueness of the Fermat point — the branch 
point position that minimizes the sum of the lengths of 
vessels at a single junction — is already well established 
(for example, see |20j). This allows us to carefully con¬ 
struct an algorithm (described in Sec. Ill C|) that reliably 
relaxes all branch point positions into the global opti¬ 
mum. 

After determining the positions of branch points for 
a given hierarchy, we compare distinct configurations to 
find the optimum network. The search through configu¬ 
rations is also a central problem in phylogenetics, where 
the goal is to construct phylogenies to identify similar 
groups of species and trace the development of genes 
through speciation. Even in the case of genes that control 
biological traits, a loosely analogous space-filling phe¬ 
nomenon emerges in the form of species filling the niches 
in the environment. With our specific goal of complete 
spatial covering of network tips, we develop strategies in 
Sec. Mi for exploring the space of hierarchies that are 
similar to those used on phylogenetic trees. 

The organization of the subsequent sections is as fol¬ 
lows. In Sec. mi we describe the basic assumptions for 
our space-filling network model, including the details of 
the local optimization of branch point locations and the 
global paths through the network. In Sec. Hill we in¬ 
troduce the specific quantitative network properties that 
we use to compare the fitness of candidate networks. We 
introduce the properties of the space of tree hierarchies 
and our implemented exploration strategies in Sec. [iVl 
In Sec. 0 we detail the results from the several layers of 
optimization that we implement, and discuss the insights 
that they offer in Sec. IVI1 


II. CONSTRUCTION OF ARTIFICIAL 
VASCULAR NETWORKS 

To better understand the connection between the local 
asymmetries of individual vessel lengths and the global 
constraint on space-filling capillaries, we optimize candi¬ 
date networks that are embedded in two spatial dimen¬ 
sions (2 -D) in silico with respect to specific optimization 
principles. We explore these optimized artificial networks 
and quantify their branching length asymmetries to com¬ 
pare with our empirical data. Our model’s simplification 
of the cardiovascular network focuses on the lengths of 
segments as defined by the straight line between adja¬ 
cent branch points. Reticulated structures occur within 
leaves to mitigate damage HEUzl and within animals as 
anastomoses (or pathologically as fistulas). However, we 
focus on the vast majority of the cardiovascular system 
that distributes resources through a hierarchical, tree-like 




3 


structure, in which no segments or subtrees rejoin down¬ 
stream to form closed loops before reaching the capillar¬ 
ies. This is sufficient for the focus of our investigation of 
the asymmetric, space-filling structure that distributes 
the resource-rich blood from the heart throughout the 
body. 

The space-filling property of the cardiovascular net¬ 
work constrains the hierarchical structure of the network 
and the positions of branch points. Here we describe 
our process for the construction of individual networks 
and the space of possible networks through the follow¬ 
ing steps: defining a distribution of space-filling service 
volumes in the space of the body, identifying all unique 
hierarchies and pairings between tips of hierarchies and 
distinct service volumes, and determining the positions 
of branch points for each hierarchy and pairing. 


A. Space-filling service volumes 

Because real systems do not organize or grow on a reg¬ 
ular (symmetric, isotropic) grid, we position service vol¬ 
umes randomly within the space they fill. Construction 
of service volumes begins by choosing a random point 
within the body volume that represents the location of a 
capillary. We then randomly choose other points (capil¬ 
lary locations) so that none he within a predefined con¬ 
stant distance from another capillary location. After de¬ 
termining a set of capillary locations that span the 2-D 
area, the entire body is partitioned into Voronoi cells fed 
by the closest capillary. In this way, each capillary be¬ 
comes associated with a specific service volume, and the 
sum of the service volumes fills the entire space (see Fig. 
[3]or[5j). 


B. Space of hierarchically distinct trees and 
pairings with service volumes 

Because multiple branching levels connect the service 
volumes to the heart, there are many possible hierarchi¬ 
cal orderings of branching junctions across these differ¬ 
ent levels. For example, there are two unique hierarchies 
when there are four service volumes: the top three con¬ 
figurations in Fig. [Id] are the same perfectly balanced 
hierarchy, while the remaining trees have the same unbal¬ 
anced hierarchy. The distinguishing feature is the pair¬ 
ings of tips in the hierarchy to specific service volumes 
(1-4). 

There are many distinct pairings of terminal tips with 
the set of fixed service volumes for each bifurcating tree. 
For four service volumes, there are a manageable 15 
unique bifurcating trees (shown in Fig. EH)- Before de¬ 
termining branch point positions for small networks, we 
exhaustively search through all possible hierarchies and 
pairings with service volumes. For large networks, the 
number of distinct hierarchies and pairings ((2n —3)!! for 


n distinct service volumes) is prohibitively large, so we 
sample the space as described in Sec. IV B1 

We do not disqualify configurations if one vessel path 
crosses with another (these would likely separate in three 
dimensions), and there is no exchange of resources or in¬ 
teraction in blood flow at such locations. Crossings are 
not observed for networks that minimize only total net¬ 
work length without a constraint on hierarchical balance, 
but they often occur for optimal configurations with a 
strong constraint on hierarchical balance. 


C. Optimization of branch point positions for a 
fixed hierarchy and pairing 

We now detail our algorithm for the optimization of 
the positions of branch points that connect a distribu¬ 
tion of service volumes to the heart through a hierarchi¬ 
cal network. Within our algorithm, the position of each 
branch point depends solely on the location of the adja¬ 
cent branch points in the network. Distant vessels affect 
each other indirectly, but not through any direct long- 
range process. Using the limited, local information given 
by the neighborhood of a branch point, each junction is 
assigned a uniquely-defined position that minimizes the 
sum of the Euclidean distances to each neighboring junc¬ 
tion. This is equivalent to the Fermat point of the tri¬ 
angle formed by the two downstream ends of the child 
vessels and the one upstream end of the parent vessel 
(see Fig. [10. The Fermat point of a triangle is a special 
case of the more general geometric median , the unique 
point that minimizes the sum of distances to an arbi¬ 
trary number of other fixed points. We follow the algo¬ 
rithm presented in [21] to avoid errors in determining the 
geometric median. Assigning branch point positions as 
geometric medians effectively minimizes the construction 
costs for the local network structure. 

We construct our networks from simple bifurcations, 
but using the Fermat point to assign branch point po¬ 
sitions can lead to coincident (degenerate) bifurcations, 
as shown in Fig. [l^ Degenerate branch points are con¬ 
solidated at the geometric median of the upstream end¬ 
point of the parent and three or more downstream end¬ 
points of the associated children segments. In this way, 
two degenerate bifurcations become a trifurcation and, 
more generally, n degenerate bifurcations become a sin¬ 
gle (n + l)-furcation. Networks that are hierarchically 
distinct in their bifurcating structure can become iden¬ 
tical networks by collapsing bifurcations. Through ex¬ 
haustive explorations (described in Sec. IV A[) . we find 
that this marginally reduces the number of possible con¬ 
figurations (see results in Fig. 112a)) . Because we have no 
a priori filter to identify which bifurcating trees are re¬ 
dundant, we must consider the entirety of the space of la¬ 
beled, rooted, bifurcating trees throughout the algorithm 
to identify sufficiently optimal configurations. With po¬ 
sitions defined in a consistent way, we can now compare 
the properties of distinct hierarchies to determine which 
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Global Optimization via Search through 
the Space of Hierarchies 
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Figure 1: (a) Schematic of the simplest bifurcating tree network, showing the heart (hollow red triangle) and two 
service volumes (filled brown circles) with labels for the lengths of each associated segment, (b) Three possible 
locations for a bifurcation junction. The rightmost configuration shows the Fermat point of A123 that minimizes 
the sum of segment lengths, (c) The two distinct bifurcations (left) collapse to a single trifurcation (center) and set 
to the geometric median of the four endpoints (right), (d) Comparing C (a measure of some length property of the 
network) for each of the 15 configurations for four service volumes shows that the bottom right configuration is 
optimal with respect to C. The relative order and position of both the tips and branch points in Fig. [id] do not 

correspond to relative positions in space. 


is the best for a particular space-filling strategy. 

III. SELECTION CRITERIA FOR BIOLOGICAL 
NETWORKS 

All characteristics of an organism that affect fitness 
and are heritable are under selection. A key question is 
which features of the vascular network are under selec¬ 
tion. Here we define specific fitness measures that are 
tied to the structure of the network configuration that 
allows us to rank candidate networks and determine the 
optimal configuration. 

A. Global length properties of space-filling 
configurations 

Here we introduce a standardized fitness measure that 
allows us to compare candidate networks for their suit¬ 
ability to transport blood and resources. Each indepen¬ 


dent measure for a network’s fitness relates to a physical 
quantity that likely guides the evolution of the cardio¬ 
vascular system toward a more efficient network. Specif¬ 
ically, the system’s cost in construction material and the 
maintained volume of the blood relates to the total net¬ 
work length — the sum of the lengths of all segments. 
Competing with this minimization of materials, the dis¬ 
sipation of the heart’s pumping power relates to the path 
lengths between each capillary and the heart. The power 
dissipated by smooth (Poiseuille) flow through a segment 
is directly proportional to the length of the segment jUj . 
In the absence of radius information, reducing the cost of 
pumping blood is equivalent to reducing the total path 
lengths that blood travels. 

We define these two fitness measures — one dealing 
with total network length and the other with individual 
path lengths — as 

L= Y, ^ (!) 

all segments i 
in network 
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H = £ £ * 

all paths p all segments i 
in network in path p 


(2) IV. GLOBAL OPTIMIZATION IN THE SPACE 
OF HIERARCHIES 


where N tips is the number of distinct service volumes, 
corresponding to the number of tips and distinct paths. 
The generalized total cost function C linearly combines 
these two measures by their respective weights Cl and 

Ch : 


C(C l ,Ch) = C l L + C h H. (3) 


This cost function connects minimization of material and 
power dissipation to study optimal networks that are 
space-filling. 

Because an increase in cost corresponds to a decrease in 
fitness, we place this approach in an evolutionary frame¬ 
work by inverting and scaling our cost function C to be 
a relative fitness function F: 


F(C Ll C H ) 


C n 


C(C L ,C H ) 


( 4 ) 


where C m i n is the most optimal network under consider¬ 
ation. By this definition, the optimal configuration has 
a fitness F = 1 and less optimal configurations have a 
fitness F < 1. 


To determine the optimal hierarchy and its connectiv¬ 
ity, we search the space of rooted, labeled, bifurcating 
trees. The positions of the branch points are fixed by the 
process in Sec. Ill Cl The globally optimal network of all 
configurations maximizes the fitness F while satisfying 
the space-filling constraint on service volumes. As an ex¬ 
ample, the optimal configuration in Fig. fTdl corresponds 
to the hierarchy in the bottom right, where the fitness 
F( 1, 0) = 1 includes only total network length L [Eq. 
©]. resulting in a Steiner tree @). 

Our exploration of configuration space has many simi¬ 
larities to phylogenetic trees, for which software is avail¬ 
able to search through the space of hierarchies mm. 
Since the available software is not tailored to our spe¬ 
cific goals of optimizing space-filling networks, we imple¬ 
ment our own algorithms. Because of the large number 
of distinct bifurcating rooted trees (that grows factori- 
ally with size), efficient search strategies generally focus 
on regions with greater fitness. We develop strategies 
to search through possible configurations and find space¬ 
filling networks that best satisfy the general biological 
constraints from Sec. urn 


B. Equal Distribution of Resources through 
Hierarchical Balance 

Because the network tends to exhibit nearly symmetric 
branching in radius and must distribute resources equally 
to each capillary in the body, the network hierarchy can¬ 
not be overly unbalanced, with one segment having many 
more tips to supply than its sibling. In accordance with 
this argument, empirical data do not show major arter¬ 
ies branching directly into capillaries. We address this 
constraint by selecting for networks with more balanced 
hierarchies. 

A hierarchy is better balanced if there are roughly 
equal numbers of tips supplied downstream by each sib¬ 
ling segment. Conversely, a hierarchy becomes more 
poorly balanced as the disparity grows between the num¬ 
ber of tips. In this sense, we define the degree that a 
hierarchy is unbalanced V as 


U = 1 — min 

all sibling 
pairs (i,j) 


J 


( 5 ) 


where N$ is the number of distinct downstream service 
volumes supplied through segment i and segment j is 
always the sibling with the most downstream tips. In 
our algorithm, we select against hierarchies for which the 
degree of unbalance U is greater than some threshold Uq. 
We eliminate configurations above this threshold before 
optimizing branch points and calculating fitness. 


A. Navigating in the space of hierarchies 


Our numerical technique guides the search by select¬ 
ing changes that increase configuration fitness. Making 
small changes in the branching structure, such as a single 
swap of two segments in the hierarchy or a regraft of one 
segment to a spatially near part of the tree, yields new 
configurations. Because the change to the hierarchy is 
small, using the positions of branch points from the pre¬ 
vious configuration saves time in optimizing the global 
positions in the new configuration. 

For local swaps in the hierarchy, we exchange a seg¬ 
ment with one child (including the associated down¬ 
stream subtrees) of the segment’s sibling (i.e., its nib- 
ling). There are 2 (n — 2) possible nibling swaps for 
n (> 2) discrete service volumes. However, nibling 
swaps do not address changes for segments that are dis¬ 
tant in the hierarchy but have small spatial separation. 
To account for these changes, we regraft single segments 
to spatially near branches of the hierarchy. We limit 
the search of spatially proximal branch points to those 
within twice the minimum service volume separation of 
each other. This restriction maintains a linear increase 
in the number of explored regrafts with the number of 
service volumes, in contrast to the factorial increase that 
would result from including all possible regrafts. 
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B. Seed for trajectories: Balanced Hierarchy 
Construction 

We accelerate the identification of near-optimal net¬ 
works by choosing an initial configuration that avoids 
many suboptimal structures (e.g. configurations with 
many repeated crossings or non-contiguous subtrees). To 
improve overall computation time, some approaches ex¬ 
plore permutations of pairing tips with service volumes 
under a constant hierarchy [26j. Fortunately, the di¬ 
mensionality of the space for each branch point posi¬ 
tion never exceeds three in our study, which allows us to 
construct a favorable configuration directly through spa¬ 
tial partitioning. Such a favorable configuration avoids 
less-fit configurations and satisfies the intuitive guidelines 
that branch points connect nearby subtrees (efficiency by 
proximity) and that sibling subtrees have similar num¬ 
bers of service volumes in accordance with symmetric 
branching in radius. 

To ensure the maximal hierarchical balance for the 
seed, we begin with a single set that contains all ter¬ 
minal service volumes. This set is then partitioned into 
two subsets of equal size (or within 1 service volume if 
the number is odd, which guarantees U < 0.5), using 
a straight line to define the boundary between the two 
sets. When appropriate, this line passes through the ge¬ 
ometric center (i.e., the unweighted average position) of 
the previous set of points and the geometric center of the 
new subset. Resulting sets further split into smaller sub¬ 
sets to yield a complete, bifurcating hierarchy. We refer 
to this process and the resulting seed as the Balanced 
Hierarchy Construction (BHC). 


C. Efficient search trajectories 

We further accelerate our search by limiting the num¬ 
ber of nearby configurations considered at each step. We 
accomplish this through a carefully guided greedy search 
through the space of hierarchies (effectively simulated an¬ 
nealing |27|, |28j at zero temperature), which often finds 
a near-optimal configuration. A greedy strategy offers 
expedited elimination of configurations that are far from 
optimal; the algorithm abandons configurations that sat¬ 
urate at a fitness lower than the current most optimal 
configuration. Our implementation allows five iterations 
of the process in Sec. ana then excludes configurations 
that fail to reduce the cost C [Eq.©] by at least 5% of 
the remaining difference from the current optimal config¬ 
uration. The algorithm with this exclusion scheme suc¬ 
cessfully identifies near-optimal configurations. 

Because the sampling process is not exhaustive, the 
search through the space of possible hierarchies is not 
guaranteed to yield a globally optimal configuration. 
However, performing a reasonably thorough search as 
we outline here and conducting several runs from the 
BHC seed (in our simulations, at least 10 runs) increase 
the likelihood of identifying a configuration that is near- 


optimal and shares many of the branching length asym¬ 
metries that an optimal configuration exhibits. With de¬ 
pendable algorithms for determining branch point posi¬ 
tions and for exploring the space of possible hierarchies, 
we can now investigate the length properties of space¬ 
filling networks under several basic space-filling strate¬ 
gies. 


V. RESULTS AND ANALYSIS 

We now present the results of optimized networks and 
of the analysis on real vascular networks, including the 
properties of the most optimal networks. To build in¬ 
tuition about the space of hierarchies, we first explore 
the space exhaustively for small networks and establish 
the distinct patterns that the two optimizations L and H 
produce. In comparing optimal configurations with ob¬ 
servations of real systems, we find better agreement by 
enforcing a constraint on the degree of unbalance U in 
the hierarchies of candidate configurations. 

A. Exhaustive search for small networks 

To become more familiar with the landscape of pos¬ 
sible configurations, we exhaustively explore the space 
of hierarchies and pairings for networks that are small 
enough to quickly yield comprehensive results for a sin¬ 
gle realization of fixed service volumes. We collapse and 
reorganize the higher-dimensional space of branch point 
swaps into a single dimension by ranking each configura¬ 
tion based on the fitness F. This reorganization involves 
a normalization of rank so that the fittest configuration 
occurs at 0 and the least fit occurs at 1. Rescaling the 
rank is necessary even for networks of the same size and 
shape because different realizations may have different 
numbers of unique configurations after consolidating de¬ 
generate bifurcations (mentioned in Sec. Ill C[) , despite 
having the same number of service volumes. Larger net¬ 
works tend to have greater range with respect to both 
costs in L and H. 

The minimum distance between each service volume 
is constant for each of the networks that constitute the 
ensemble of realizations for the curves in Fig. [2a] Each 
curve represents the average fitness (relative to the fittest 
configuration for the particular set of capillary positions) 
over an ensemble of networks with a fixed number of tips 
and constant total area. In generating the ensemble, we 
exclude those that arise with a different number of tips 
than what is desired until we accumulate 1000 configu¬ 
rations of the target size. Across curves, the total area 
increases to produce networks with more service volumes 
more frequently. 

One might expect a large set of similarly fit, near- 
optimal networks, which would be represented by a 
plateau near the optimum. However, the sharp descent 
away from the optimal configuration in Fig. I2al indicates 
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(a) (b) 

Figure 2: (Color) Ranked landscapes for unique consolidated configurations from exhaustive exploration of 
bifurcating trees of fixed area and fixed Nu ps (a) Exhaustive landscapes for fitness based only on total network 
length L [F( 1, 0)]. (b) Exhaustive landscapes for fitness based total network length L and average path length H 
[F( 1, 9)]. We choose the weights (Cl, Ch ) = (1,9) so that the contribution from H is not dominated by the 

contribution from L. 


that there are few configurations that are near-optimal. 
From an evolutionary perspective, this implies that the 
vascular networks of organisms are under strong selec¬ 
tion. Furthermore, the slope near the optimum becomes 
steeper as more service volumes are introduced, so that 
the best configurations become more distinct from other 
possibilities as the number of service volumes grows. 
Considering the very large number of service volumes 
in real organisms, this again indicates that real vascular 
networks are under strong selection pressures for space¬ 
filling and efficiency. 

Optimal networks that have no constraint on hierarchi¬ 
cal balance fall into two general classifications depending 
on the relative weights of total network length L and av¬ 
erage path length H in the fitness measure F. As shown 
in the simple examples of Fig. [3j network fitness mea¬ 
sures that are weighted to minimize L yield bifurcating 
trees, while measures that are weighted to minimize H 
yield hubs. Bifurcating trees better correspond to real 
networks, suggesting that total network length L plays 
a larger role than average path length H. Since a sin¬ 
gle hub is not observed (and not expected from material 
costs) in real systems, we do not consider configurations 
that ignore total network length L. However, optimizing 
only for L leads to meandering, bifurcating paths, which 
become shorter and more direct by including both costs 
(L and H). Furthermore, additional global information is 
necessary to directly minimize path lengths than the local 


environment that we consider in Sec. Ill Cl — specifically, 
the context of the entire path. This means that our anal¬ 
ysis is best suited for optimality that always includes a 
significant contribution from total network length L and 
a weaker contribution from average path length H. 


B. Trajectories for sampling larger networks 

With better intuition about the space of hierarchies 
from small networks, we now explore the space for larger 
networks with more service volumes. The branching 
properties of larger networks give more applicable results 
to connect particular space-filling strategies with the ob¬ 
servations of real cardiovascular systems. We first sum¬ 
marize the properties of optimized networks without any 
constraint on hierarchical balance (Uq = 1). 

Because the search through the space of hierarchies is 
not exhaustive for large networks, we cannot show ranked 
landscapes averaged over ensembles with different service 
volume positions as we did for small networks. Instead, 
we show landscapes from a single realization of service 
volume positions that come from an ensemble of trajec¬ 
tories that start with the BHC configuration and end at 
a local optimum (Fig. Illal in App. [Bj . The greedy al¬ 
gorithm samples fewer less-fit configurations, yielding a 
shallower slope near the optimum than the exhaustively 
explored landscapes in Fig. [2a] Since the starting point 


















Global Optimization for 
Total Network Length L 

[J’(l.O)] 



Global Optimization for 
Average Path Length H 

[f(0,l)] 



Figure 3: (Color) Two classes of networks: (a) The 
optimal configuration that minimizes total network 
length [L in Eq. (pQ)] of the 15 possible trees 
(corresponding to Fig. EU consists only of bifurcations, 
(b) The optimal configuration that minimizes average 
path length between each service volume and the heart 
[H from Eq. ©] of the 15 possible trees consists of a 
single hub. The regions of varying background color 
define the Voronoi cells corresponding to individual 
service volumes. 



Figure 4: (Color) Average fitness landscapes for total 
network length L [F( 1, 0)] over 100 trajectories for a 
single network of each size. 


of the search (the BHC configuration) is already favor¬ 
able, we expect that the worst-ranked configuration of 
the partial search is already very near optimal. 

Searches through the space of hierarchies and the prop¬ 
erties of optimal configurations do not vary with differ¬ 
ent convex body shapes. Fig. [5] shows example opti¬ 
mized networks for a maximally symmetric body shape 
(see App. [B] for other shapes). The general trends of 
long, meandering paths for solely minimizing total net¬ 
work length L and of more direct paths when including 
H are consistent across both isotropic, circular areas and 
elongated, rectangular areas. 

To characterize branching features of these large con¬ 
figurations, we quantify the asymmetric branching at¬ 
tributes with the two ratios 

Al = (6) 

^C2 

\R= V f (7) 

choosing £ Cl < £ C2 for the lengths of child 1 and 2 and 
r c i < r C2 for the radii (shown in Fig. [laj)- Note that 
perfect symmetry corresponds to Al = A r = 1 and 
smaller values of Al and A r correspond to more asym¬ 
metric branching. Distributions for the branching asym¬ 
metry ratio in length Al for various sizes and shapes 
are shown in Fig. [6l Branch points for which one of 
the child segments does not exist because of degener¬ 
acy with a service volume center do not contribute to 


the distribution for Al- There is little change in the 
features of the distribution of Al across different sizes 
of networks with Uo = 1. This trend persists for both 
isotropic and anisotropic enclosing shapes (as Fig. [6] in 
App. 151 shows). A summary for the cross-generational 
length ratio is given in App. [Cl 


Many branch points in these networks coincide with 
a service volume, predicting large trunks that feed cap¬ 
illaries directly. Similar results appear in the study of 
flow through a dynamic, adaptive network [29). How¬ 
ever, such a trend does not agree with the empirical data. 
Although there is asymmetry in adjacent segments at 
branch points and a lack of strict balance in the hierar¬ 
chy along different paths, we observe that large arteries 
do not branch directly to capillaries and arrive at the 
same expectation from the dynamics of blood flow. The 
major qualitative distinction between the BHC and the 
optimized configurations with Uq = 1 is that the BHC 
is a network with a balanced hierarchy. Upon inspection 
of empirical data in Sec. Ivcl we find that the branch¬ 
ing length asymmetries for the BHC configuration (given 
in App. EJ motivate an additional constraint on hier¬ 
archical balance during the search through the space of 
hierarchies. 
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BHC Seed 


Total Network Length L [F(l, 0)] 


Total Network Length L and Mean 
Path Length H [F(l, 9)] 



Figure 5: (Color) Optimal configurations for two fitness measures. We choose the weights (Cl, Cr ) = (1,9) so 
that the contribution from H is not dominated by the contribution from L. 



Figure 6: (Color) Distributions of Al for several sizes of 
circular areas, averaged over 100 realizations of service 
volume distributions and optimized solely for total 
network length L [F( 1, 0)]. Local length asymmetry 
between siblings skews toward symmetry for optimal 
networks with no hierarchical balance. 


C. Comparison of optimized networks with 
empirical data 

The results in Sec. IV Bl show that optimization for 
total network length or average path length with no con¬ 
straint on hierarchical balance leads to distributions of 


asymmetry in sibling vessel length that skew toward sym¬ 
metry (Al ~ 1). We now present the analysis of Al 
that characterizes the local length asymmetries at branch 
points for real and optimized networks. From this anal¬ 
ysis, we explore how limiting the degree of unbalance U 
in an optimal artificial network yields asymmetries that 
better match biological networks. 


1. Asymmetric vessel length distributions of real networks 

We analyze MRI images of the human head and torso 
as well as micro tomography images from wild-type 
mouse lung. Both data sets break from strict symmetry. 
As shown in Fig. [7J the network-wide distribution for 
A r is skewed toward symmetry (A r ~ 1), while the dis¬ 
tribution for Al is more uniform, representing a greater 
contribution from very asymmetric branching (Al < 1). 
These results are representative of general features for 
length distributions in real biological networks. The fact 
that the optimized networks in Sec. IV Bl do not exhibit 
a similar distribution for Al signals that important bio¬ 
logical factors are missing. Because of the skew toward 
symmetry in sibling segment radii, we limit the hierar¬ 
chical unbalance of optimized networks in Sec. IV C 21 


2. Degree of balance necessary to match biological networks 

Imposing a constraint on hierarchical balance leads to 
configurations that reflect more realistic asymmetry in 
branching lengths. Hierarchical balance, which equalizes 
the number of service volumes that each sibling segment 
supplies, is related to the blood flow that is required 
to deliver resources and effectively limits the asymme- 
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Sibling Radius Ratio (A r) 


Sibling Length Ratio (A l) 




(a) (b) 

Figure 7: (Color) Observed radius and length branching asymmetry ratios [Eqs. (|6]) and (0, respectively] in mouse 
lung and human torsos, (a) Radii ratios are skewed toward symmetry (A r ~ 1), although they are not always 
perfectly symmetric, (b) Length ratios are not skewed toward symmetry (many ratios have A l < 1), contrary to 

symmetric models. 


try of sibling radii. In Fig. [8] we show results for sev¬ 
eral thresholds for the constraint on hierarchical balance. 
Decreasing the threshold Uo yields more realistic distri¬ 
butions for Aas shown in Fig. [9] Because these net¬ 
works are embedded in 2-D, decreasing Uo can also result 
in more crossings between segments at different levels. 
By comparing Figs. I9aland l9bl we see that the constraint 
on hierarchical balance leads to similar results indepen¬ 
dent of the weight of average path length H in config¬ 
uration fitness. Instead of contributing significantly to 
fitness, H is effectively optimized through hierarchical 
balance. 


While enforcing hierarchical balance leads to more re¬ 
alistic branching and length asymmetry distributions, it 
is not necessary to have a maximally balanced hierarchy. 
In FigfTOl we show that lowering the threshold Uo reduces 
network fitness. As the constraint on hierarchical balance 
Uo decreases, the average fitness also decreases. Since 
the distribution of Xl is approximately uniform around 
Uo ~ 0.7 and below, the best value for the hierarchical 
constraint is Uo ~ 0.7 because this yields the fittest net¬ 
works that have uniform distributions for length asym¬ 
metry. 


VI. DISCUSSION 

With our determination of branch point positions and 
exploration of distinct hierarchical configurations, we can 
remark on several consequences that follow from the gen¬ 
eral properties of optimized networks. Organizing the 
lengths between branch points to fill 2- or 3-D space 
with capillaries inevitably leads to asymmetries and un¬ 
balanced networks [30]. Strictly symmetric and balanced 
networks are either inefficient in materials or not space¬ 
filling. For example, in the H-tree all children branch or¬ 
thogonally from the parent, resulting in inefficient paths. 
Other networks with more efficient paths lead to capil¬ 
laries that are equidistant from the source, which could 
cover the surface of a sphere but not fill its volume. For 
the optimal, space-filling networks that we explore, we 
impose a constraint that pushes the network toward hi¬ 
erarchically balanced branching structures but does not 
require maximum balance. One can imagine other inter¬ 
esting metrics for hierarchical balance, but we concen¬ 
trate on how a maximum degree of unbalance Uo affects 
the structure of the network. This guarantees a min¬ 
imum level of balance in the hierarchy but still allows 
freedom in the search for optimal networks, as well as 
nonuniformity in the hierarchical balance. 

We construct a seed configuration that builds a net¬ 
work to ensure maximal hierarchical balance while main¬ 
taining efficient contiguity of subtrees. Configurations 
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Uo = 0.6 

Figure 8: (Color) Optimal configurations for several 
constraints on hierarchical balance and two 
optimization weights for fitness F. The BHC seed is the 
same as in Fig. [5] We choose the weights 
(Cl, Ch ) = (1,9) so that the contribution from H is 
not dominated by the contribution from L. 


that tend to be hierarchically balanced, such as the BHC 
configuration (where the constraint is implicit in the 
construction algorithm) or optimized configurations that 
limit unbalance, do not show a strong skew toward sym¬ 
metric branching in lengths. This hierarchical balance 
may result from gradual, incremental growth as an indi¬ 
vidual organism matures and ages. Nearby vessels grow 
to supply resources to new tissue, resulting in contiguous 
subtrees and favoring routes that reduce path lengths and 
avoid a single, meandering artery that branches directly 
to capillaries. 

Other computational models approach the growth and 
optimization of space-filling networks in different ways. 
Although there are many algorithms to generate struc¬ 
ture that do not intentionally optimize network archi¬ 
tecture or space-filling properties, near-optimal configu¬ 
rations may emerge spontaneously from certain simple 


Because of the many different factors and interactions 
that influence the structure of the cardiovascular system, 
our basic model can be expanded in many directions. 
Radius information can be incorporated into optimized 
networks by requiring flow to be uniform in all terminal 
service volumes. By including radius information, blood 
flow as well as more appropriate structural and energetic 
costs can lead to revised optimization principles, which 
require the calculation of the weighted Fermat point (e.g., 
see [go]) and has been explored previously in a limited, lo¬ 
cal context [H,[33, [38|. Note that lowering the threshold 
Uo tends to increase the minimum number of branching 
levels between the heart and capillaries. Less drastic hier¬ 
archical unbalance implies that the ratios of parent-child 
radii f3 = r c /r p should be near 1 (symmetric branching). 
This translates the global, topological property into a 
local branching quantity. 

We do not expect that increasing the dimensionality of 
our networks to 3-D would change the qualitative results 
for branching asymmetry in length (Al) with hierarchical 
balance (specifically Uo)- However, the numerical loca¬ 
tion for an optimal trade-off between fitness and balance 
may shift. Studies of large vessels (near the heart) show 
these vessels to be planar 39 }), but the planarity cannot 
always hold across the entire network if tips must fill a 
3-D space. Still, in the absence of obstacles, all opti¬ 
mization conditions enforce planarity in 3-D for branch 
points in their local context. Introducing regions where 
the network is prohibited (e.g. through bones, organs or 
from self-avoidance) constrains the Fermat point to the 
surface of a sphere or some other shape [40, [4Jj . 

While the topological change of allowing loops intro¬ 
duces many complications to the properties of flow and 
hierarchical labels 142], such a modification can be ben¬ 
eficial in understanding reticulated vascular structures. 
Loops are especially important when considering network 



Uo = 1.0 


Uo = 0.9 


Total Network Length L 

m, o)] 


rules. Examples of such pattern formation processes and 
Total Network Length lassociated algorithmic rules include models for both an- 
and Mean Path Lengthgiogenesis (18l . l3l| . a s well as vasculogenesis (in terms 
H [F( 1, 9)] Q f chemotaxic [32L l33j|. mechanical substratum Q, and 

cellular Potts models [35, l36j|). However, these models 
do not adequately address our focus on branching length 
asymmetries for efficient, hierarchical, space-filling net¬ 
works. Specifically, the pattern formation model for an¬ 
giogenesis does not incorporate consistent space-filling 
service volumes, only space-filling arterial structure. The 
arterial structure fills some regions so that they are de¬ 
void of capillaries, while multiple tips converge to the 
same location elsewhere. The models for vasculogenesis 
do not optimize the development of a hierarchical branch¬ 
ing network. However, dynamic vascular remodeling [29] 
can form structures both with and without closed loops 
while maintaining a uniform distribution of capillaries, 
although the optimal structures also suffer from large ar¬ 
teries branching directly to capillaries. We extend these 
models to understand the asymmetric lengths of adja¬ 
cent segments in vascular networks and how these relate 
to space-filling service volumes. 
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A l for Total Network 
Length L [F( 1, 0)] 


A l for Total Network 
Length L and Average Path 
Lpnp-th H \F(1. Qh 
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Figure 9: (Color) Distributions for several thresholds of hierarchical unbalance Uo- (a) Distributions of Xl optimized 
solely for total network length L [F(1,0)]. (b) Distributions for A l optimized for total network length L and average 
path lengths H [F(l,9)]. All plots are averaged over 200 realizations of service volume distributions. 



Figure 10: (Color) Average fitnesses F for ensembles of 
100 optimized configurations. Configurations are more 
fit if a greater hierarchical unbalance Uo is allowed. 


robustness (i.e., resilience to damage) within organs and 
leaves [l6|, Q3] ) or pathological growth in tumors Ml ii- 
These types of network properties can be included in fu¬ 
ture models. 


Locally, the position of a branching junction minimizes 
the sum of vessel lengths in our model. Globally, we im¬ 
pose a threshold on the minimum hierarchical balance, 
which reduces the differential blood flow into sibling seg¬ 
ments. Although real vascular networks consist almost 
entirely of bifurcations (although there is rapid, asym¬ 
metric branching from the aorta to capillaries through 
coronary arteries), the iterative approach described in 
Sec. Ill Cl can lead to low numbers of bifurcating junc¬ 
tions for some candidate networks. 

Limiting the degree of unbalance in the hierarchy does 
not continue to shift the distribution of A^ away from 
symmetry (A l ~ 1) below Uq « 0.7, which suggests 
that there is an appropriate trade-off between the hierar¬ 
chical balance threshold Uo and configuration fitness F 
that does not require perfect symmetry for an efficient 
network structure. The increased cost of the network 
in Fig. [10] is similar for both curves, implying that the 
increase mostly comes from total network length L. 

The large number of distinct bifurcating hierarchies ne¬ 
cessitates that we carefully choose and execute the algo¬ 
rithms for searching the space of possible configurations. 
Consequently, we construct a favorable starting point and 
concentrate computational resources on regions that are 
most likely to contain optimal configurations. Using the 
numerical implementations in sections ITT1 and HVl we iden¬ 
tify optimal networks and study the length properties of 
individual segments within the context of a network with 
space-filling terminal service volumes. 

Our results have many implications for how vascular 
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networks fill space efficiently. We exhaustively explore 
fitness landscapes for small networks and carefully guide 
the sampling of the space of hierarchies for large net¬ 
works in order to determine near-optimal configurations. 
Our results show that strict hierarchical balance is not 
optimal for the architecture of cardiovascular networks. 
Furthermore, there is a trade-off between hierarchical 
balance (which is related to symmetric branching in ra¬ 
dius at the local level) and the distribution for branch¬ 
ing in lengths that shows the connection between the 
space-filling and efficiency requirements of the network. 
By incorporating radius and flow information, as well as 
growth patterns that incorporate obstacles and loops, we 
can continue to build on present models to better under¬ 
stand vascular architecture and gain insights for its ef¬ 
fects on resource delivery, metabolic scaling, aging, and 
repair after damage. 


Appendix A: Similarity measure to compare 
hierarchical groupings between configurations 

While collapsing the landscape of measures to a single 
dimension informs us about the typical distribution of 
configurations, it retains no information about the rela¬ 
tion of the hierarchies between different trees. To address 
this issue, we define a measure of similarity to compare 
how two hierarchies group the same set of tips. This 
measure is normalized such that similar hierarchies and 
groupings of service volumes have a similarity score near 
1 , while hierarchies that group service volumes in very 
different ways have a similarity score near 0. To meet 
these guidelines, we perform a simple count of the num¬ 
ber of identical subtree groupings between two hierar¬ 
chies and normalize by the maximum possible number 
that could be shared if the trees were identical. In accor¬ 
dance with these properties, define the similarity S(A, B ) 
between two configurations A and B as 

or A TD\ — __ 

1 5 j “ma x{a(A,A),a(H,H)} 

<W') = £ £ (fmCn T 4cm) 

subtree m subtree n 

in network X in network Y 

where I s is the indicator function (1 if statement s is true 
and 0 otherwise) and subtree refers to the set of tips in 
that particular subtree. 

Configurations that have a worse fitness measure are 
less similar to the optimal configuration, as shown in Fig. 
nn However, note that similarity S(Ti,T opt ) is not a 
monotonic function when rank i is defined by the con¬ 
figuration’s measure. For example, consider hierarchy A, 
which may be very similar to a hierarchy B, which it¬ 
self is very similar to C. Then it is possible that A and 


C are less similar to each other than each is to B, yet 
both are ranked higher than B with respect to a par¬ 
ticular measure. This also means that optimal config¬ 
urations are not always a single swap or regraft away 
from all near-optimal configurations, i.e. local minima 
are possible. The average similarity in Fig. Illbl does not 
approach 1 , indicating that the subtree grouping of ser¬ 
vice volumes can be very different between networks that 
are nearly optimal (F « F opt ). Some of the stratification 
into distinct levels of similarity is apparent in Fig. Illal 
for the single trajectories (the “Single BHC Seed” and 
“Best Random Seed”). 

Appendix B: Network size and shape 

In Fig. I12al we show the distribution of the number of 
distinct hierarchies after consolidating degenerate bifur¬ 
cations for these ensembles of fixed numbers of service 
volumes. The variance in the number of unique config¬ 
urations increases with network size, but a dominating 
contribution that increases the number of configurations 
comes from adding a service volume. Not surprisingly, 
the average total number of service volumes (. N t i ps ) scales 
linearly with the total area (see Fig. I12bj) . In Fig. [14] 
we show distributions for branching length asymmetry 
A l. All distributions for optimized networks that have 
no constraint on the hierarchy (i.e. Uq = 1) exhibit a 
strong skew toward symmetry. 

Appendix C: Parent-child length ratio 7 

In Fig. I15al we show the network-wide distribution for 
the ratio of child-to-parent lengths 7 = i c !i p for child c 
with parent p. Although there is the tendency that 7 < 
1 , some child segments have a relatively shorter parent. 
Although slight, an increased threshold for Uo shifts more 
child segments to be shorter than their associated parent 
(see Fig. rnvhii . Independent of the threshold, the nonzero 
variance of this distribution shows that 7 is not constant 
throughout the network. 
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Figure 11: (Color) Partial search trajectories in the the space of hierarchies, (a) Ranked fitness landscapes for the 
partially explored hierarchical space during optimization for total network length L [F( 1, 0)]. (b) Corresponding 
similarity measures. All quantities are calculated with respect to the most optimal configuration found over all runs 
The solid line shows the average fitness over 100 trajectories from the BHC for a system with 36 service volumes. 
Individual paths are for a network of 48 service volumes. One path begins with the BHC configuration and the 
other begins with a random seed (the best performing of 5 random seeds). 



(a) 



(b) 


Figure 12: (Color) Network properties with increased size, (a) The number of unique configurations N con f ig for each 
fixed number of service volumes is narrowly distributed relative to the increased number of configurations from 
introducing an additional service volume, (b) The average number of service volumes increases directly proportional 
to the total area.The solid grey line represents the standard deviation about the mean. 
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Figure 13: (Color) Optimal configurations for several shapes and two fitness measures. Both isotropic (circular) and 
anisotropic (rectangular) areas exhibit similar space-filling strategies for near-optimal configurations. We choose the 
weights (Cl, Ch) = (1,9) so that the contribution from H is not dominated by the contribution from L. 
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Figure 15: (Color) Distributions for the ratio of child-to-parent length 7 . (a) Distributions of 7 in mouse lung and 
human head and torso, (b) Distributions of 7 for several thresholds of hierarchical unbalance [To, optimized solely 

for total network length L [F(1,0)]. 
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